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Abstract. An orbit tracing technique was used to generate current sheets for three magne- 
totail models. Groups of ions were followed to calculate the resulting cross-tail current. 
Several groups then were combined to produce a current sheet. The goal is a model in 
which the ions and associated electrons carry the electric current distribution needed to 
generate the magnetic field B in which ion orbits were traced. The region -20 /?£ < x < -14 
Rg in geocentric solar magnetospheric coordinates was studied. Emphasis was placed on 
identifying the categories of ion orbits which contribute most to the cross-tail current and 
on gaining physical insight into the manner by which the ions carry the observed current 
distribution. Ions that were trapped near z = 0, ions that magnetically mirrored throughout 
the current sheet, and ions that mirrored near the Earth all were needed. The current sheet 
structure was determined primarily by ion magnetization currents. Electrons of the 
observed energies carried relatively little cross-tail current in these quiet time current 
sheets. Distribution functions were generated and integrated to evaluate fluid parameters. 
An earlier model in which B depended only on z produced a consistent current sheet, but it 
did not provide a realistic representation of the Earth’s middle magnetotail. In the present 
study, B changed substantially in the x and z directions but only weakly in the y direction 
within our region of interest. Plasmas with three characteristic particle energies were used 
with each of the magnetic field models. A plasma was found for each model in which the 
density, average energy, cross-tail current, and bulk flow velocity agreed well with satel- 
lite observations. 


1. Introduction 

Current sheets are present in the Earth’s magnetotail, magneto- 
pause and bow shock, and in many other space plasmas. The adia- 
baticity parameter is defined as tc 2 = R c / p zo , where R c is the 
magnetic field line radius of curvature and p zo is the gyroradius of 
an average current-carrying particle at z = 0, the center of the cur- 
rent sheet [ Buchner and Zelenyi , 1989]. The criterion tc < 1 for the 
principal current-carrying particles may be used to define a thin 
current sheet. Such particles do not follow simple guiding center 
orbits. The principal distinction between the sheet-like structures 
studied here and a ring current is that the magnetic field in a current 
sheet must be generated by particles within the sheet. In contrast, 
the Earth’s internal dipole dominates the magnetic field in the radi- 
ation belts. It therefore is necessary to consider local plasma cur- 
rents to model the current sheet, while it may be possible to neglect 
local currents in the radiation belt. The goal of the present work is 
to investigate the physical structure of the thin current sheet in the 
Earth’s magnetotail. 

Consistent Orbit Tracing (COT) 

Several methods have been used to construct models of the 
magnetotail. Each method uses a different set of basic assumptions 
and therefore provides different and often complementary insights 
into the physical structure of a current sheet. Simulations and some 

Copyright 1996 by the American Geophysical Union. 

Paper number 96JA01945. 

0148-0227/96/96JA-01945$09.00 


orbit tracing methods involve injecting particles into a spatial 
region with fixed boundaries. Knowledge therefore is needed about 
the source of particles and the boundary conditions. Initial mag- 
netic and electric field configurations within the region of interest 
are used to start a simulation, but detailed knowledge of these 
fields is only needed at the boundaries. Electric and magnetic fields 
are calculated self-consistently within the region of interest as the 
system evolves in time. 

In contrast, the COT method requires that the basic magnetic 
and electric fields are known throughout region of interest. How- 
ever, very little is assumed about particle sources or particle bound- 
ary conditions. Our selection of the COT technique was based on 
the assumption that the global magnetic field models give a good 
description of the magnetotail but that less is known about particle 
sources. One major limitation is that only steady state conditions 
can be studied. 

The COT technique is equivalent to solving the Vlasov equa- 
tion, so involves the motion of a non interacting particle in a field 
generated by all other particles. Wave-particle interactions and 
instabilities are not considered in this treatment. The results pre- 
sented here are confined to the region -20 R E < x < -14 R % , 0 < Izl 
< 2 R e in geocentric solar magnetospheric coordinates. This region 
of interest was selected because AMPTE/IRM and ISEE 1 data 
were available. The major task in the COT procedure is to find a 
distribution of ions and electrons that will as nearly as possible 
generate the preselected magnetic field. In some respects the 
method is similar to a BGK solution [Bernstein et al. t 1957]. The 
BGK method is based on the assumption that a particle distribution 
can be found that will, through Poisson’s equation, generate a wide 
variety of preselected electric field structures. The COT method is 
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based on the assumption that a particle distribution can be found 
that will, through Ampere’s law, generate a variety of preselected 
magnetic field structures. If a realistic tail magnetic field model has 
been picked, then we assume that some distribution of ions and 
electrons with energies and densities similar to those measured in 
the tail can be found to produce this field. 

The first step was to follow many ion orbits in the original mag- 
netic and electric fields. The use of an error limit of 10 -10 for tra- 
jectory tracing assures the inclusion of all nonguiding center 
effects. This detailed inclusion of orbital information is one of the 
strengths of the COT technique. Each group in the present study 
consisted of 1000 ions. Earthward drift was produced by a uniform 
dawn to dusk electric field. All current sheet ions were traced as 
they came in from tail ward of the region of interest and drifted 
earthward completely through this region. No physical boundaries 
exist in the system except for the very rare particles that hit the 
Earth. 

The cases described here started with monoeneigetic ions. The 
angular characteristics of a biMaxwellian were used to randomly 
select an initial phase and pitch angle for each ion in the group. Ion 
energies change as particles drift in the electric field, but ions in 
each group remain in a fairly narrow energy band. Monoenergetic 
ions were used for this particular study because each such group 
tends to be dominated by particles that follow one of the orbit cate- 
gories that are described below. A principal goal of the study was 
to find which orbit types are needed to produce a typical magneto- 
tail current sheet. When Maxwellian or other broadly distributed 
ion groups are used, each group contains a mixture of ions follow- 
ing many different types of orbits. It therefore is harder to gain 
physical insight into the kinds of orbits that are needed to generate 
each model magnetic field. 

As described later, a number of assumptions were introduced in 
order to add electron currents. The ion plus electron groups were 
then combined using a linear least squares fitting routine to find 
particles which will generate the preselected magnetic field. 

One-Dimensional Results 

The COT technique was used by Kaufmann and Lu [1993] to 
construct an approximately consistent one-dimensional (1-D) 
model magnetotail. Although orbits were traced in three dimen- 
sions, the model is referred to as a 1-D magnetotail because the 
modified Harris current sheet magnetic field was used [ Harris , 
1962]. The modified Harris model consists of a constant B zo and a 
component B x (z) = B^ 0 tanh(z/L) that varies only in the z 
direction. 

Even though nearly consistent 1-D current sheets were found, 
the model did not provide a reasonable approximation to the mag- 
netotail in our region of interest. The problems with the modified 
Harris model were that it required a density of several-keV ions 
that was higher than is observed, and it required streaming or 
anisotropies at the edge of the current sheet that were much larger 
than those observed. These are inherent problems with any consis- 
tent 1-D model of the middle magnetotail. An unrealistically high 
density is needed because there is no net cross-tail drift of trapped 
particles in a 1-D magnetic field [Stern and Palmadesso , 1975]. We 
were unable to generate modified Harris current sheets with typical 
magnetotail field strengths using ions and electrons with the ener- 
gies observed unless the density was 1 cm -3 or larger. The average 
observed density in the -20 R E <x < -M R E region of interest is 
approximately 0.3 cm” 3 [Baumjohann et al. f 1989; Huang and 
Frank , 1994b]. 


It is well known that streaming or a strong anisotropy, near the 
firehose instability limit, is needed at the edge of any 1-D current 
sheet to produce force balance [Rich et al. y 1972; Cowley , 1978]. 
The difficulty is that 1 -D models have no x dependence, so require 
anisotropies much larger than those observed to balance the j x B 
force or equivalently, in a consistent model, the tension force of the 
sharply curved magnetic field lines. For these reasons, we con- 
cluded that no 1-D current sheet can produce a reasonable model of 
the middle magnetotail. 

The present paper shows that the above problems are not 
present in a magnetic field model which has both x and z depen- 
dence. Section 2 describes the magnetic field models used, section 

3 describes the ion orbits that exist in the model fields, and section 

4 presents the basic COT results. Those fluid parameters that can 
be most directly compared to satellite measurements are discussed 
in section 5. Finally, section 6 summarizes the conclusions. Future 
papers will present the model distribution functions, analyze more 
complex fluid characteristics such as force balance and energy 
flow, and consider the question of uniqueness. 

2. Tail Model 

A simple magnetic field model was developed that provides a 
good representation of the portion of the tail being studied. A sim- 
plified model was needed to make accurate tracing of many ion 
orbits practical. The present study used a sum of fields from an 
Earth centered dipole, a ring current [Tsyganenko and Usmanov , 
1982], an equilibrium or Harris-like tail field [Birn et al. y 1975; 
Zwingmann , 1983], and a uniform B zo . Three-dimensional versions 
of the dipole and ring current fields were used so the mirroring of 
particles at low altitudes is as realistic as possible. The equilibrium 
tail field is 2-D, with B varying in the x and z directions. 

The equilibrium tail field module is dominant in the region of 
interest. This module represents magnetic fields that are generated 
by approximately isotropic particle distributions. The specific form 
used was the one designated ZW by Karimabadi et al . [1990]. The 
dipole and uniform B zo are vacuum fields, so require no electric 
currents in the magnetotail. The ring current is largely carried by 
particles earthward of the region of interest. 

The dipole field was fixed for all versions of the model. Each of 
the other modules had parameters that were adjusted to generate a 
magnetotail with desired properties. Adjustments first were made 
so that the k parameter and the total integrated cross-tail current 
would be similar to those at midnight in the Kp = 4 version of the 
much more complicated T89 [ Tsyganenko , 1989] magnetosphere. 
An empirical magnetotail model tends to have a thicker current 
sheet than is present in the actual magnetotail. Waves and flapping 
of the tail broaden the apparent observed thickness obtained by 
binning satellite data. McComas et ai [1986], Sanny et al. [1994], 
and Zhou et al. [1995], using the ISEE 1 and 2 satellite pair, con- 
cluded that the typical instantaneous scale size of the principal cur- 
rent sheet is of the order of 1 R E in our region of interest. We 
selected the Kp = 4 version of T89 to adjust k because this version 
is thinner and has a larger equatorial B z than the smaller Kp ver- 
sions. 

Figure 1 shows field lines in the noon-midnight meridian plane 
for the T89 model and the three simpler models used in the present 
study. Within the region of interest, characteristic scale lengths in 
the z direction for the cross-tail current j y are about 0.3-0.4, 0.7, 
and 1.0 R e in the models labeled thin, standard, and thick, respec- 
tively. The corresponding T89 scale length is about 1.2 R E . In each 
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Figure 1. Magnetic field lines in the noon-midnight meridian plane for the (a) T89 with Kp = 4, (b) thin, (c) standard, 
and (d) thick magnetotail models. 


case the current sheet thickness decreases by about 10% as one 
moves in from x = -20 to -14 R E . 

The extra uniform B z0 term was included in our model so the 
magnetic field at z = 0 could be adjusted. For the standard model, 
B z at the equator dropped from 5 nT at x = -14 R E to below 3 nT at 
x = -20 R e . To get the same k, the equatorial B z dropped from 6 to 
about 3.5 nT in the thin model and from 4.5 to slightly above 2 nT 
in the thick model. This parameter is important in determining k 
and the type of orbit that is dominant in any region of the tail. 

Some uncertainty remains concerning the most appropriate 
equatorial B z to select. Huang and Frank [1994a] defined the neu- 
tral sheet as the region with \B X \ < 5 nT in ISEE 1 data from 1978 
and 1979. They found the average B z to be approximately 7 nT in 
the region -22 R E <x< -10 R E . Fairfield [1986] used many years 
of data from the IMP 6, 7, and 8 satellites taken within ±3 R E of the 
estimated location of the neutral sheet. Data were divided into one 
group for -20 R E < x < -10 R E and a second group that extended 
from x = -20 R e to the apogee of each satellite. This study found 
B z as large as 7 nT to be uncommon in the lower altitude range, and 
extremely rare beyond 20 R E . Rostoker and Skone [1993] used 
consistent methods to select only clear neutral sheet crossings in 
ISEE 1 data for -22 R E < x < -10 R E and in IMP 8 data for -38 R E 
<x< -22 R e . Their Figure 3b shows a discontinuous drop in B zo at 
x = -22 /?£, which they suggested may partly be attributable to 
observations at different phases of the solar cycle. A corresponding 
average B z has not been determined for the 1985 and 1986 
AMPTE/IRM data sets we have used. Nakamura et al. [1991] 
selected a group of neutral sheet crossings by requiring |B| to drop 
below 5 nT. Since the observations with which we will make initial 
distribution function comparisons were selected from this group of 
neutral sheet crossings, most have B z in the range that was used in 
the simple models. However, Nakamura et al. [1991 ] noted that the 
estimated uncertainty of the measured B z is 1 nT while other com- 
ponents of B have uncertainties of only 0.1 nT in AMPTE/IRM 
data. This is because the satellite spin axis, along which magnetic 
field calibration is most difficult, lies close to the solar magneto- 
spheric z direction. 

3. Orbit I^pes 

Ion orbits in the magnetotail usually are classified according to 
the particle’s dynamical characteristics, for example, a chaotic, a 


resonant, or a guiding center orbit [Chen, 1992]. Particle dynamics 
are strongly dependent upon k. However, the cross-tail current pat- 
tern carried by ions in the k range studied here was found to be 
more dependent upon the particle’s mirror point location than upon 
the dynamic orbital characteristic. Since a primary goal was to 
study the structure of a consistent current sheet, this paper catego- 
rizes orbits into three groups according to mirror points. We first 
define z Q = mv/ [qB x (z Q ) ] , which is the point at which the radius 
of curvature in the y-z plane of a particle with mass m, charge q, 
and velocity v is equal to the particle’s distance from the equatorial 
plane. The first distinctive category consists of meandering and 
other trapped particles, which return to the equatorial plane without 
moving beyond Izl = 2 z 0 . Particles that magnetically mirror within 
the principal current sheet at 2 z 0 < Izl < 2 R e comprise the second 
category. These are sometimes called cucumber orbits [Buchner 
and Zelenyi , 1989]. Finally, the third category involves particles 
that mirror closer to the Earth. These are sometimes called Speiser 
orbits. The following section will show that some contribution 
from each of the above three orbit categories was needed to create 
a consistent current sheet that realistically approximated the 
Earth’s magnetotail. 

Sample Orbits 

The sample orbits in this section show the physical manner in 
which particles carry cross-tail current in the standard tail model 
(Figure lc). Magnetization currents are dominant for particles that 
mirror inside the current sheet [Bird and Beard , 1972; Kaufmann 
and Lu , 1993]. The x-z and x-y projections of the proton trajectory 
in Figures 2a and 2b show all three orbit categories. This proton 
had 5 keV of energy and was located at (x, y, z) = (-18.5, 0, 0.5) R E 
when tracing started. The ion mirrored once within the northern 
hemisphere of the current sheet as it drifted earthward from the 
starting point, became briefly trapped at Izl < 2 z Q = 0.5 R E near x = 
-16 to -18 /?£, and then mirrored three times in the southern hemi- 
sphere. The first of these three southern hemisphere mirror points 
was within the principal current sheet and the other two were closer 
to the Earth. The orbit then was traced backward in time from the 
starting point, showing that this proton also had been trapped 
between x = -29 and x = -20 R E . This alternation between periods 
of being trapped near z = 0 and periods of magnetically mirroring 
at various points is typical of chaotic orbits [Ashour-Abdalla et al 
1995]. 
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Figure 2. (a), (b) Two projections of one proton orbit in the stan- 
dard magnetotail model. Tracing started with a 5-keV proton at (*, 
y y z ) = (-18.5, 0, 0.5) R e . (c), (d) Tracing started with a 5-keV pro- 
ton at (-16.5, 0, 1.5) R E . The arrows point to the starting locations. 


the details of a trajectory and of deviations from guiding center 
motion. Figure 2 also shows the most important differences 
between orbits in 1-D and either 2-D or 3-D models. Note that dur- 
ing the last two orbit segments in Figure 2b, when the mirror point 
was below jc = — 10 R E> the particle followed almost the same field 
line as it moved down to and back from a mirror point. The same 
feature is evident in Figure 2d whenever the ion mirrored below 10 
R e . Protons mirroring near the Earth with energies of about 5 keV 
in this magnetic field appear to simply bounce off the current sheet 
[Speiser, 1965] when they reach the equator (Figure 2a). Mirror 
points on the orbits immediately before and after a given current 
sheet interaction are nearly the same if the particle’s k parameter is 
close to a resonant value [Buchner and Zelenyi , 1989; Chen> 1992]. 
The mirror points of chaotic particles, with tc well removed from 
the resonant values, generally are different after each current sheet 
interaction. However, in this paper all resonant and chaotic parti- 
cles that mirror near the Earth are placed into the same category 
because the cross-tail current distribution in the region of interest is 
similar for these two dynamically different trajectory classes. 
Whenever the mirror point is near Earth, the total cross-tail dis- 
placement achieved in one bounce cycle is approximately equal to 
2 gyroradii based on B zo . Note that different scales are used for the 
x, y, and z axes in Figure 2, so that circular motion near z = 0 
appears to be stretched in the y direction. The motion of particles 
that mirror at low altitudes in 2-D and 3-D models is similar to the 
motion of untrapped particles in a 1-D model. In all these cases, j^ 
is concentrated in a sheet near z = 0 that is much thinner than is 
required to consistently generate the standard model magnetotail. 

Figure 2b shows that this ion drifted part way back in the nega- 
tive y direction when it mirrored within the current sheet, near z = 
-1.5 R e> x = -10 R e and near z = 0.7 R E , x = -14 R e . In both cases 
there was a substantial net cross-tail drift during a full bounce 
cycle. Similar bounce cycles are seen in Figure 2d. The fact that all 
trapped and mirroring ions undergo a net cross-tail drift during 
each bounce cycle is one of the most important differences 
between either 2-D or 3-D and 1-D models. In the 1-D models, 
where both the x and y components of the canonical momentum are 
conserved, all mirroring particles drift equal distances in the +y 
direction near z = 0 and in the -y direction as they mirror. Self-con- 
sistent 2-D and 3-D models have lower total plasma densities than 
1-D models because particles that mirror within the current sheet 
carry a net cross-tail current in the x dependent models. 

Particles of species s that obey the guiding center equations in a 
gyrotropic plasma carry a current given by 


Since the -20 R E < x < -14 /? £ , 0 < Izl < 2 R E region was being 
investigated, orbit tracing was stopped on the earthward side when 
the ion crossed z = 0 at a point at least 2 gyroradii earthward of x = 
-14 R e . Similarly, tracing backward in time was stopped when the 
ion crossed x = -20 R E at least 2 gyroradii beyond Izl = 2 R E . An 
equatorial crossing tailward of the field line that connects to x = 
-20 /?£, Izl = 2 R e plus 2 gyroradii also can be used to stop back 
tracing. It is evident from Figure 2 that no specific particle bound- 
ary conditions were imposed at the edges of the region of interest, 
and that nothing unusual happened at the starting position. Starting 
points were selected to get groups of ions that were dominated by a 
particular trajectory type. For example, starting the ions with an 
isotropic distribution at z = 0 gives a group that begins with many 
trapped ions, while no ions that start at large Izl are initially trapped. 

Cross-Tail Drift and Current 

Figures 2b and 2d show the physical manner by which the 
bounce averaged cross- tail drift becomes relatively independent of 
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This current involves both drifts of the guiding centers and 
a divergence-free magnetization drift v m . The small polariza- 
tion drift has been neglected. The drift \ g is a combination of 
gradient, curvature, and ExB drifts. Curvature guiding cen- 
ter drift is associated with the factor in (1). The curvature 
guiding center drift is large and positive near z = 0. The first 
P s± factor in (1) represents temperature and density gradient 
magnetization currents, and the second P sl factor in (1) is the 
curvature magnetization current or orbit crowding effect. The 
magnetic field gradient produces both a drift of guiding cen- 
ters in the +y direction near z = 0 and a magnetization current 
in the -y direction. No magnetic field gradient term appears in 
(1) because the associated guiding center and magnetiza- 
tion \ m terms cancel. 
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Figure 2b shows that there is a small cross-tail drift per current 
sheet interaction during the periods when the ion was trapped at Izl 
< 0.5 /?£, beyond x = —20 R E . This helps to explain why the 
bounce-averaged drift velocity depends relatively weakly on orbit 
type. Trapped particles make many current sheet interactions, but 
have a net cross-tail displacement per interaction that is much less 
than 2 gyroradii. Particles mirroring near the Earth have few cur- 
rent sheet interactions, but have a net displacement of nearly 2 full 
gyroradii during each interaction. 

Figures 2c and 2d show another important property of magneto- 
tail current carriers. This trajectory is characterized by a series of 
events in which the ion mirrors and then interacts briefly with the 
current sheet. The cross-tail motion tends to increase modestly dur- 
ing each current sheet interaction as the ion moves into the inner 
tail because the ion is gaining energy. However, it is more impor- 
tant that the earthward drift speed is decreasing at low altitudes. 
The ion moves earthward by as much as 10 R E during the time 
required to mirror when its equatorial crossing point is near x = 
-30 to -40 /?£. Earthward motion is only about 1 R E per mirror 
cycle near x - -10 R E . This decrease in earthward drift is the pri- 
mary reason that the displacement between adjacent equatorial 
crossing points is mostly in the x direction in the distant tail and 
mostly in the y direction in the near-Earth tail. The result is an 
Alfv6n layer or wall [Ashour-Abdalla et al . , 1992]. A group of par- 
ticles from a deep- tail source with a given average magnetic 
moment do not penetrate to lower field lines. 

Figures 3a, 3b, and 3c show x-z, x-y , and y-z plots of a proton 
that started at (-15.5, 0, 1.0) R E with 15 keV of energy. This high 
energy proton passed through both the N = 1 and N = 2 resonances. 



Figure 3. Three projections of one proton orbit in the standard 
magnetotail model. Tracing started with a 15-keV proton at (-15.5, 
0, 1.0 )R e . 



gies change as the particles drift. 


Resonances are characterized by symmetric orbits that involve tit- 
tle net change in magnetic moment. A particle’s magnetic moment 
changes dramatically near the halfway point between resonances. 
A resonant ion crosses the z = 0 plane N+l times during each cur- 
rent sheet interaction [ Buchner and 2telenyi, 1989], An individual 
particle in the tail is accelerated as it drifts earthward and 
duskward. Acceleration tends to increase and therefore to 
reduce K for the particle. However, k tends to increase when the 
particle reaches the more dipolar field lines with large R c in the 
inner tail. Figure 4 shows that these two competing effects lead to a 
minimum k near x = — 17 R E in our models for all starting eneigies 
used. The k parameter of the particle in Figure 3 decreased from a 
little more than 0.5 at the x = -32 R E crossing, which lies near the 
N = 1 resonance, to about 0.3 at -15 /? £ , which is near the N = 2 
resonance. Note that the ion crossed z = 0 twice during each of the 
two most distant neutral sheet interactions before it became briefly 
trapped, and then crossed z = 0 three times during each of the four 
most earthward interactions. The magnetic moment and the mirror 
magnetic field changed substantially during the crossings at -25 R E 
< x < ~20 R e when k was near 0.4. This behavior is typical of 
highly chaotic motion halfway between resonances. 
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Figures 5a, 5b, and 5c show x-z, x-y, and y-z trajectory projec- 
tions for a proton that started with 15 keV of energy at (-14.5, 0, 
0). Although it is unusual to last this long, the orbit shows that it is 
possible for a proton to remain trapped in a very simple pattern as 
it drifts earthward all the way from x = -30 R E to x = -15 R E . The 
K parameter changed substantially in this interval. This orbit is 
referred to as a figure eight pattern because of its appearance in the 
y-z projection. By following one cycle of the figure eight pattern in 
Figure 5c, it is particularly easy to see that this ion always moves in 
the +y direction at Izl >0.3 R E and in the —y direction at Izl < z 0 . 
Although not as obvious, all groups of trapped ions carry a similar 
current distribution. 

4. Generation of Model Current Sheets 

This section presents sample results from nine runs using the 
COT technique. The thin, standard, and thick magnetotail models 
were each used for three runs. For each magnetotail model, one 
plasma sheet was generated using monoenergetic 1.5 keV, one 
using 5 keV, and one using 15 keV starting proton energies. Ions 
are energized as they drift in the uniform Ey — 0.3 mV/m or 1 .9 
kV/J? £ electric field. However, the ions in an initially monoener- 
getic group remain concentrated in a relatively narrow energy band 
at any one location, so that a single orbit type tends to dominate. 

Distribution Functions 

The calculations used twenty 0. \-R E thick z boxes in the range 0 
< Izl < 2 R e . The model is symmetric about z = 0, so orbital infor- 



Figure 5. Three projections of one proton orbit in the standard 
magnetotail model. Tracing started with a 15-keV proton at (-14.5, 
0, 0 )R e . 


mation at equal distances above and below z = 0 were combined in 
these 20 boxes. Six 1 -R E wide x boxes in the range -20 R E <x< 
-14 R e were used at each z. Since there is very little y dependence 
in this region of interest, only one y box extending toy = ±10 R E 
was used. As a result, distribution functions were generated for a 
total of 120 spatial boxes. We refer to this set of runs as two dimen- 
sional because of the lack of significant y dependence even though 
all calculations and binning are fully three dimensional. Although 
no useful information would be gained by keeping multiple y boxes 
for the present runs, several y boxes are retained when either y 
dependence is added to the equilibrium tail module or when distri- 
bution functions are needed near or sunward of x = -10 R E . 

Two methods were used to calculate the ion number density 
«■ (x, z) and current density j, ( x , z) . The point at which each 
box was entered was determined accurately during orbit tracing. 
New orbit points were introduced at each box edge crossing point 
so that all orbit segments remained within one of the 120 spatial 
boxes. Use of the efficient Bulirsch-Stoer tracing algorithm made 
step sizes so large that linear interpolation was not adequate for the 
curved trajectories. The new techniques described below removed 
the anomalies at certain pitch and phase angles that were noted by 
Kaufmann and Lu [1993]. The most serious problem in our earlier 
study involved the generation of distribution functions by counting 
particles only when they crossed a plane surface. With large step 
sizes, the distribution functions generated by this method had 
anomalously low values when the velocity normal to the detector 
was small. 

The first technique involved keeping track of the total time 
spent and the total distance that all particles in a group traveled in 
the x, y, and z directions while they were inside each box. Density 
is proportional to this time and j ( - (x, z) is proportional to the dis- 
tance. Eastwood [1972] and Kaufmann and Lu [1993] described 
the use of this method. The resulting j; (x, z) includes all ion guid- 
ing center drift, magnetization, and nonguiding center currents. 

For the second technique, velocity space was divided into 30 x 
30 x 30 boxes in the v,, v y , and v z directions. The result is a 
27,000-point velocity distribution function for each of the 120 spa- 
tial boxes. The distribution functions were generated by adding the 
time spent when taking each step along a trajectory into two of the 
27,000 x 120 velocity-spatial boxes, with half the time for each 
step placed into the initial and half into the final velocity box for 
this step. Integrating the distribution function gave essentially the 
same n ( (x, z) and j ( - (x, z) as was obtained by the first method. 
Adding in the electron current, described below, gave the total 
j (x, z) in each box. 

Sample Case 

A set of 30 groups of 1000 ions were traced for each of the nine 
cases noted above. The angular distributions all were selected so 
that /(r, v) would have been isotropic if the ions had started in a 
uniform B. An actual /( r, v) may be anisotropic even in the start- 
ing box because ions can return to this box with different pitch 
angles. For any one of the nine runs, the 30 groups of ions all 
started with the same energy, which was 1 .5 keV, 5 keV, or 1 5 ke V. 
The only difference between the 30 groups in any one run was the 
location at which orbit tracing began. The starting points were z = 
0, 0.5, 1.0, and 1 .5 R e at the centers of each of the six x boxes, plus 
at z = 2.0 on field lines that crossed the equator at the centers of 
each x box. 

The three columns in Figure 6 show the cross-tail components 
of the ion, electron, and total ion -plus -electron group currents for 
the thick-5 keV case. Figures 6a through 6e show currents from 
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Selected Groups for Thick Sheet, 5 keV 


Ion j y Electron j y Group j y 



Figure 6* (a) - (e) Cross tail currents carried by groups of 1000 ions, their accompanying electrons, and by the total 
ion-plus-electron groups. The thick magnetotail model was used, (f) Weighted sum of Figures 6a to 6e. The solid 
curve in the group j y column is the current needed to produce a consistent magnetotail. 


five groups of 1000 particles, and Figure 6f is a weighted sum of A number of assumptions were made to calculate the current 
Figures 6a through 6e. For example, the ion jy panel in Figure 6a carried by Maxwellian electrons. Most of these assumptions have 

shows the current carried by a group of 1000 ions that were started already been described because they were used in the 1-D calcula- 

at x = -15.5 R E , z = 0. The smooth solid curve in the group jy col- tions [Kaufmann and Lu , 1993]. A polynomial fit was made to both 

umn of Figure 6f shows the cross-tail current that is needed to gen- n / (x, z) and the average ion energy for each group traced. The 

erate the model magnetic field. electron number density n e (jc, z) was set equal to the polynomial 

The horizontal axis of each panel in this figure is divided by tick fit to n t (jc, z) so the resulting plasma would be approximately 

marks into six blocks associated with the six different x boxes. The neutral. Since n e would be constant along field lines if electrons 

block at the left side of each panel, labeled -19.5, includes all were isotropic at the equator, a first-order parallel electric field was 

information for the range -20 R E <x< -19 R E . The block at the added to the zeroth-order uniform E y . The Boltzmann relation 

right side of each panel, labeled -14.5, corresponds to -15 R E < x < 

-14 R e . Within each x block are 20 points representing increasing z n e ( x , z) = n e (x 0 , z 0 ) exp [~q e $/T e ] (2) 

values in the range 0 < Izl < 2 R E , The smooth curve in the group jy 

column of Figure 6f shows that the j y (x, z) which is needed to was used to calculate the electric potential difference 0 

generate B is largest at z * 0 in each x box, and drops monotoni- between each (x, z) box and a reference point (x 0 , z 0 ) on the 

cally as z increases to 2 R E . This panel also shows that larger cur- same field line. The average electron energy or 3/2 T e was 

rents are needed in x boxes closer to Earth. The 5-keV row in taken to be 1/7 of a polynomial fit to the ion group average 

Figure 7 shows the same information as is in the group jy panel in energy [Bawnjoharm et a/., 1989]. Since electrons are much 

Figure 6f. The six x boxes are plotted separately and labeled more less energetic than ions, electrons were influenced much more 

completely in Figure 7. strongly by the parallel electric field than were ions. Equation 
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Figure 7. Solid curves are the cross-tail currents need to produce the thick model current sheet. Each column corre- 
sponds to a l-#£ thick x box. Dashed curves are the least squares fits obtained by adding currents carried by five to 
nine groups of 1000 ions and electrons. There are 20 points along the horizontal axis in each panel, corresponding to 
the 20 z boxes. The three rows of panels were produced by starting protons with 1.5, 5, or 15 keV of energy, as indi- 
cated. 


(2) only gives the distribution of along a field line, so only 
determines the parallel electric field. Since E u is distributed 
differently along adjacent field lines, some E ± in the x-z plane 
must also be present. We used a variational approach to deter- 
mine this E ± . It was assumed that E ± will develop so that the 
energy density e 0 E 2 is minimized. This assumption com- 
pletes the definition of E. The total potential drop between z = 
0 and z = 2 R E was no more than 200 V. This potential differ- 
ence is similar to the values obtained in our previous 1-D 
model. Some groups of ion orbits were traced again in the 
magnetic field with this added £ t) . The small potential drop 
had very little effect on the ion j y (x, z) distribution. 

The K for each group of electrons is 1 1 times the k for the asso- 
ciated proton group, so all electron motion is well approximated by 
the guiding center assumptions in the models studied here. The 
electron cross-field current was calculated from (1) based on the 
polynomial fits to n e (x, z) and an isotropic P e (x, z) = n e T e . 
This procedure yielded the electron j y (x, z) for each group, as 
shown in the middle column of Figure 6. The addition of electron 
and ion currents gave the total j y (x, z) for each group, as shown 
in the right column of Figure 6. 

The above procedure, which assumes isotropy and the presence 
of £j| , is not the only one that could be used to include electrons. 
Instead, it is possible to produce charge neutrality along field lines 
with £|| = 0 by adjusting the angular distribution of electrons at 
the equator. Substantial electron anisotropies involving T, j > T L 
have been observed in the central plasma sheet [Hada et ai t 1981; 
Paterson et a/., 1995], though most observations have reported 
approximate isotropy at z = 0. In the guiding center approximation, 


a biMaxwellian distribution at the equator yields a density varia- 
tion along a field line of 

" <e| - (« y » i [ -«* /r “ 1 <3) 

where n 0 and B 0 are the particle density and IB I at the equator. 
When £|| = 0 or <|> = 0 everywhere, the density varies by a fac- 
tor of 7*1, /T ± as one moves from the equator to the iono- 
sphere. 

When the particle density is known all along a field line, it is 
possible to calculate the angular distribution needed for charge 
neutrality with £|| = 0 throughout the region where the guiding 
center equations are valid. However, we retained density informa- 
tion only in the previously defined region of interest, so do not 
have densities all the way down to the ionosphere. One model in 
which the ion density changes rapidly along a field line will be 
shown later. An electron distribution with very abrupt pitch angle 
variations would be needed to maintain charge neutrality in this 
case if E u = 0. The method used in the present analysis is easy to 
implement and results in reasonable distribution functions for all 
cases studied. Since electron cross- tail currents are small in sheets 
with the thicknesses used here, this simple method should be ade- 
quate. 

Current Distributions 

It is clear that none of the individual groups shown in Figure 
6a-6e carried a current distribution that is even remotely similar to 
the current needed for self-consistency. This observation was valid 
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for each of the 30 groups of particles traced in each of the nine 
cases studied. However, a weighted sum of several groups often 
can produce a reasonable approximation to the necessary 
j y (x, z) . The final current density obtained by linearly combining 
the five groups in Figure 6a-6e is shown by the dashed line in the 
group jy panel of Figure 6f and in the 5-keV row of Figure 7. 

Current distributions produced by the individual ion groups in 
the ion j y column of Figure 6 show the typical patterns that were 
described qualitatively in the discussion of individual orbits. For 
example, many figure eight and other trapped orbits appear in the 
group of particles shown in Figure 6a, which started at z = 0 at * = 
-15.5 /?£. This is the second of the six x blocks from the right. The 
resulting ion current shows the typical trapped particle signature 
involving negative j y at z = 0 and positive j y beyond z Q% which is 
near or below 0.3 R E for 5 keV protons in the region of interest. 
The above current pattern is dominant in the starting and adjacent x 
boxes (three rightmost blocks in the panel) because many particles 
remain trapped for only a few current sheet interactions. Similarly, 
the ion j y panel in Figure 6c shows the current carried by ions that 
started at z = 0, x = -17.5 R E . The characteristic trapped particle 
j y (x, z) pattern here appears in the more distant x boxes (blocks 
on the left side). 

There is a net drift of all trapped ions in the positive y direction, 
as noted previously. The net positive current can be seen in Figures 
6a and 6c because the area under the positive j y spikes is a little 
larger than the area under the negative spikes. However, the struc- 
ture of j y (x, z) clearly is dominated by the sharp positive and 
negative peaks, which are associated with magnetization currents, 
rather than by a small relatively uniform positive j y at Izl < 0.5 R E 
that would be predicted by (1) for an isotropic distribution. 

Ions in Figure 6e were started at z = 2 R E on a field line that 
crosses the equator at x = -19.5 R E . These ions all have initial mir- 
ror points outside the current sheet. The ions carry cross-tail cur- 
rent that is strongly confined near z = 0 when they first bounce off 
the current sheet, in x boxes near the left side of the ion j y panel in 
Figure 6e. The current is less well confined during succeeding 
bounces at lower altitudes because some ions mirrored closer to z = 
0 or became trapped. In addition, ions with low altitude mirror 
points usually have drifted several Earth radii earthward before 
they return to z = 0. This j y (x, z) structure is typical of any group 
of particles that is selected to be generally isotropic far from z = 0, 
regardless of the particle’s k or dynamical characteristics. The 
principal difference between particles on chaotic and resonant 
orbits is that resonant particles would maintain low mirror points 
for several bounces, so the pattern of narrow regions with strong 
positive j y would be more extended in the x direction. A chaotic 
particle group would not continue to be dominated by particles 
with low mirror points after one or two current sheet interactions. 

Figures 6b and 6d show currents carried by ions that started 
with an isotropic distribution at z = 0.5 %, with x = -15.5 and 
-17.5 R E respectively. These groups are dominated by ions which 
mirror somewhere in the outer portion of the current sheet, beyond 
z = 0.5 R e . The associated j y patterns tend to be complex, but are 
reproduced by other randomly selected groups of 1000 ions, even 
down to small details. Trapped ions with mirror points throughout 
the current sheet are needed to broaden the current carried by 
5-keV protons so they can produce a sheet that is thicker than the 
region associated with meandering orbits. 

As noted above, none of the individual ion groups has a j y struc- 
ture that is similar to the goal, as shown by the solid curves in the 
group j y panel of Figure 6f and in the 5-keV row of Figure 7. For 
example, the low mirror point group (Figure 6e) produces much 


too thin a current sheet. However, if some of the trapped particle 
group currents (Figures 6a and 6c) are added, it will decrease the 
positive current very close to z = 0 and add some positive current at 
larger izl. This is the procedure used in the COT technique, as 
described previously [Kaujmann and Lu, 1993]. 

A linear combination using all 30 particle groups that were 
traced for each of the nine cases studied would yield the best least 
squares fit to the j (x, z) that is needed to generate a consistent 
magnetic field. However, a principal goal of this study was to see 
which orbit types were essential for the production of a consistent 
current sheet. For this reason, only those groups that improved the 
fit at a 95% confidence level were retained. Positive weighting 
coefficients were required for each group. Of the 30 groups gener- 
ated for each case, no more than nine and no fewer than five groups 
were selected after imposing these requirements. The thick 5-keV 
case shown in Figure 6 was used as an illustration, even though the 
final fit was one of the poorest, because it was the one with the few- 
est groups in this final fit. 

Figures 7-9 show how closely the model j (x, z) compares to 
the desired j (x, z) that is needed for self-consistency for each of 
the nine cases studied. This is the most critical step in the COT pro- 
cedure, and the validity of the technique depends upon our ability 
to find consistent solutions which also agree with observations for 
nearly all preselected current sheet structures. The existence of a 
permanent current sheet in the magnetotail suggests that it is possi- 
ble to find some such distribution of particles, though not necessar- 
ily with groups of particles which all have the same starting energy. 
The jaggedness of the nine least squares fits is a consequence of the 
summation of only five to nine nearly monoenergetic particle 
groups. Smoother results are obtained in other runs using Max- 
wellian distributions and distributed injection points. 

The 1 5-keV results in Figures 7-9 show good agreement with 
the required j y for all three magnetic field models. However, aver- 
age ion energies are observed to be closer to 5 keV in the typical 
quiet magnetotail. The thin and standard magnetic field model 
results using 5-keV starting energies yielded j y that is in reasonable 
agreement with that needed for consistency (Figures 8 and 9). 
However, we did not find 5-keV ion groups that carried j y over a 
sufficiently broad range in z to match the currents needed in the 
thick model at large Ixl. Problems in finding j (x, z) that have the 
desired z dependence also are evident in the 1. 5-keV results, partic- 
ularly beyond x = - 1 7 R E . 

Part of the fitting problem may be associated with our use of the 
equilibrium tail module as the dominant contribution to the mag- 
netic field in our region of interest. The module was selected 
because it is known that this field can be generated by an isotropic 
distribution function. However, the distribution functions from 
which the entire set of equilibrium tail modules were derived con- 
tain unrealistic properties. One property that is in disagreement 
with observations is the uniform cross-tail drift velocity at all x and 
z that is assumed when deriving the equilibrium modules [Harris, 
1962; Schindler , 1972; Kan , 1973]. In contrast, observations and 
our COT results have the largest near z = 0. Burkhart and Chen 
[1993] pointed out the unrealistic thickening of the distant tail cur- 
rent sheet in this set of models. The complete fields, which have 
been produced by selecting parameters and adding several modules 
to fit T89 (Figure 1), do not show a dramatic thickening. Neverthe- 
less, the difficulty in finding ions which will generate a thick 
enough current sheet in the more distant portion of our tail models 
may be associated with this unrealistic property of the equilibrium 
tail module. A final unrealistic assumption used when developing 
the equilibrium tail module was that a frame exists in which E = 0 




z, r e 

Figure 9. Similar to Figure 7, except using the thin magnetotail model. 
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everywhere. No such coordinate system exists in a 2-D or 3-D tail 
if E y is uniform in the Earth frame. 

A different type of problem is caused by our use of monoener- 
getic particle groups. It is hard to find any completely monoener- 
getic group that carries significant j y near Izl = z 0 . Reasonable fits 
were found in most cases because there is a substantial spread of 
energies in the final COT distribution functions, which are 
weighted sums of five to nine groups. Particles were all started 
with a given energy, but each group is started at a different loca- 
tion. The energy of each particle changes as it drifts to a desired 
observation point. 

One important result is that, in all nine cases studied, at least 
one group dominated by each of the three orbit categories that 
carry distinctive j y ( x , z) distributions had to be retained. The 
COT technique therefore was unable to generate a current sheet 
with the desired thickness using generally isotropic ions with the 
observed energies unless particles mirroring near the Earth, parti- 
cles trapped near z = 0, and particles mirroring throughout the prin- 
cipal current sheet were present. Some 1-D models have been 
produced using only untrapped particles that are injected at large 
Izl, but the generally isotropic 2-D and 3-D models appear to also 
require trapped particles and particles with cucumber orbits. 

5. Fluid Parameters 

A full ion velocity distribution function was saved for each of 
the 120 spatial boxes. Various moments of f { (r, v) were evaluated 
to calculate the ion density, bulk flow velocity, and other parame- 
ters. This section discusses the parameters that are easiest to com- 
pare with published satellite measurements. It should be 
emphasized that the least squares fitting routine that determined 
which groups to keep used only j (x, z) to set weighting factors. 
The other parameters that were preselected for each run were the 
magnetic and electric fields, the ion starting energies, an isotropic 
starting pitch angle distribution, and the 1/7 ratio of electron to ion 
energies. No other adjustments were made to any of the parameters 
discussed in this section in an attempt to obtain agreement with 
observations. 

Density and Parallel Electric Fields 

The solid lines in Figure 10a are contours of constant n i ( x , z) 
and the dashed curves are magnetic field lines. Baumjohann et al. 
[1989] binned AMPTE/IRM data according to satellite location 
and the AE index. The two spatial bins that are most relevant to our 
model of the current sheet were together referred to as the inner 
CPS. The region closest to the neutral sheet was defined by B < 
7.5 nT, where B 2 xy = B 2 X + B y , and will be referred to as the inner 
current sheet. The portion of the inner CPS farther from z — 0 was 
defined by 7.5 < B^ < 15 nT, and will be referred to as the outer 
current sheet. The AMPTE/IRM satellite did not make measure- 
ments at the statistically averaged location of the neutral sheet 
[Fairfield, 1980] in our region of interest during the 1986 magneto- 
tail data collection periods. However, many neutral sheet crossings 
were observed when the neutral sheet was displaced from its aver- 
age location. The entire plasma sheet is known to undergo substan- 
tial motion in the z direction. No obvious dependence of plasma 
parameters upon the neutral sheet displacement from its average 
location was noted in a brief comparison with 1985 data, which 
had better orbital coverage. The inner CPS data show an average 
density of about 0.3 cm - -* with little z dependence. 

Kistler et al. [1993] sorted AMPTE/IRM data taken at radial 
distances beyond 15 R E according to substorm phase and to dis- 




Figure 10. (a) Contour plot of lines of constant density in 
ions/cm for the standard 5-keV model. Dashed curves are mag- 
netic field lines, (b) Contour plot of the average ion energy in keV 
for the standard 5-keV model. 


tance from the neutral sheet. This study used a more complex 
method to estimate distance from the neutral sheet. The plasma p, 
or ratio of particle to magnetic field pressures was determined for 
each data point. The highest p is expected near z = 0, and the low- 
est p in the lobes. Points were binned according to p into 3 groups, 
labeled inner, middle, and outer plasma sheet. With respect to sort- 
ing by substorm phase, it is most reasonable to compare our quiet 
time model to data taken before substorm onset. Densities were 
about 0.4 cm -3 in the inner and middle third, and about 0.2 cm -3 in 
the outer third of the plasma sheet (L. Kistler, private communica- 
tion, 1995). 

Huang and Frank [1994b] sorted ISEE 1 data according to dis- 
tance from the Gosling et al. [1986] average neutral sheet location. 
The satellite orbit provided good coverage of this neutral sheet 
region. The average observed density was approximately 0.25 
cm“ 3 during both quiet and disturbed times, with little consistent z 
dependence. 

The above measurements involved statistical averages over a 
number of orbits. An absolute distance scale was not determined in 
these statistical treatments. For example, it usually was not known 
when the satellite was exactly 0.5 R E , 1.0 R E , 1.5 R £y or any other 
specific distance from the center of the neutral sheet. Only a few 
analyses have been earned out that determined the instantaneous 
structure of the current sheet during quiet conditions in our region 
of interest. These studies measured the current sheet structure 
using both ISEE 1 and 2 spacecraft. Of these, we could find only 
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two that also described plasma density and temperature variations 
within the current sheet. McComas et al. [1986] studied three 
crossings that were associated with the passage of a strong inter- 
planetary shock front. It was concluded that the plasma density, 
temperature, and pressure changes were consistent with approxi- 
mate plasma plus magnetic field pressure balance. However, tem- 
poral changes in the current sheet structure and orientation were so 
large during this orbit that it was not possible to determine a con- 
sistent density or temperature profile. Zhou et al. [1995] found two 
periods during which a number of current sheet crossings took 
place with the spacecraft well positioned to measure density and 
temperature gradients. They concluded that both density and tem- 
perature decreased by a factor of two as the spacecraft moved l .5 
current sheet thickness scale lengths from z = 0. The average half 
thicknesses on this orbit were 2 R E when the IMF pointed north- 
ward and 0.5 R E when the IMF pointed southward. 

The COT standard 5-keV model number densities nearest the 
neutral sheet (Figure 10a) are close to or slightly below the 
observed values. The dropoff of n (*, z) with increasing z is simi- 
lar to that measured by Zhou et al. [1995] but more rapid than the 
dropoff suggested by the statistical studies. The model density 
depends on the starting ion energies used. Densities for the stan- 
dard model were 3 to 4 times larger when 1.5-keV protons were 
used instead of 5-keV protons, and about half as large as in Figure 
10a when 1 5-keV protons were used. The density in the thick 
model with 5-keV protons was similar to that in Figure 10a, while 
densities near z = 0 were about 50% larger in the thin 5-keV model. 
Our earlier self-consistent 1-D model used 3-keV protons [Kauf- 
mann and Lu % 1993] and required particle densities as high as the 
2-D model with 1.5-keV protons. 

Figure 10a shows that the model ion density is nearly constant 
along magnetic field lines near z = 0, and that the density increases 
slowly as one moves earthward along a field line beyond z - 1.5 
R e . However, the density increases rapidly as one moves earth- 
ward in the 1 R E < z < 1.5 R E region. As described previously, 
charge neutrality was maintained in the model consisting of isotro- 
pic electrons and nonguiding center ions by introducing a parallel 
electric field. Very abrupt changes would be required in the equato- 
rial pitch angle distribution of electrons on these field lines if the 
equally rapid drop in electron density that is needed to maintain 
charge neutrality were to be produced by mirror effects with E u = 
0. The density changes seen along field lines at z < 1 R E and at z > 
1.5 R e in Figure 10a were typical of the changes seen throughout 
the region of interest for the other eight cases studied. No other 
case showed as large a field-aligned density drop as that seen at 1 
R e <z< 1.5 R E in Figure 10a. 

Thermal Energies 

Figure 10b shows the average ion energy for the standard 5-keV 
model. Baumjohann et al. [1989] observed an average inner 
plasma sheet temperature of 4 keV in the region of interest, with 
slightly lower values during quiet times. Huang and Frank [1994b] 
found average temperatures of 3-4 keV during quiet times, and up 
to 7 keV during disturbed times. The temperature was higher in the 
inner than in the outer plasma sheet. L. Kistler (private communi- 
cation, 1995) found temperatures dropping from 4 keV in the inner 
plasma sheet to 2 keV in the outer plasma sheet before substorms. 
Zhou et al. [1995] found approximately a factor of two drop in 
temperature when moving 1 .5 thickness scale lengths from z = 0. A 
temperature of 4 keV corresponds to an average ion energy of 3/2 
7=6 keV. Figure 10b shows average energies rather than tempera- 
tures because monoenergetic groups were used. 


The calculated average energy in the region of interest dropped 
from about 7 keV to about 3.5 keV, which is close to the observed 
temperature variation. Since the starting ion energy is one of the 
preselected parameters, the average energy in the model can easily 
be changed. As noted above, lower energies will increase the 
required density, so a small decrease in the starting energy to below 
5 keV would slightly improve the agreement with both the mea- 
sured n (jc, z) and 7 (x, z) . 

The observed plasma p or particle/field energy density ratio 
[Baumjohann et al ., 1989] dropped from approximately 20 in the 
neutral sheet region to 3 in the outer current sheet, and to 0.3 in the 
outer plasma sheet. Results from all nine of the cases studied here 
showed decreases from p = 20 - 1 00 at z = 0 to (3 = 0. 1 - 0.5 at z = 2 
r e . The detailed dropoff depended weakly on the particle energy 
and magnetic field models used. All cases therefore showed rea- 
sonable agreement with observations. 

Bulk Flow 

Finally, Figure 1 1 shows the model ion bulk flow velocity 

V ia (r) = (1/n) Jv a /, (r, y)d\ (4) 

Detailed examination of /; ( r, v) shows that two features con- 
tribute to this bulk flow. The principal cause of V iy is an asym- 
metry in the nearly monoenergetic spherical phase space shell 
of ions in our model. Ions going to the east and west have the 
same energies, but there are more going westward. This 
behavior is associated with the large westward displacement 
of ions, especially near z = 0, that is evident in Figures 2 and 
3. This net positive produces most of the cross-tail current, 
with eastward electron drift causing the rest. It is interesting to 
note that the model V iy carried by a realistic particle distribu- 
tion decreases with increasing Izl even though the dominant 
equilibrium tail field module was developed using a uniform 
V*. 

The nature of V u shown in Figure 11 is substantially different. 
This component is primarily field aligned except very near z = 0, 
where V ix is mostly E x B drift. The model distribution function 
away from z = 0 shows that the entire spherical shell experiences a 
net shift in the positive V* direction, so that ions moving earthward 
are slightly more energetic than are those moving tailward. A less 
important angular asymmetry also is sometimes seen in the oppo- 
site direction, with a larger number of ions moving tailward than 
earthward. This small angular asymmetry reduces the net earth- 
ward Vfr , particularly near z = 1 R E in Figure 1 1 . A decrease in V ^ 
somewhere near z = 1 R E was seen in most of the nine cases stud- 
ied. A positive V# is commonly seen in our region of interest in 
both AMPTE/IRM [Baumjohann et al ., 1989] and ISEE 1 [Huang 
andFranK 1994b] data. 

In a 1-D model, the uniform E y = 0.3 mV/m field is transformed 
to zero in the deHoffman Teller reference frame which moves 
earthward at a speed E y /B zo . The z component of B is uniform in a 
1-D model. There is no reference frame in which this uniform E y is 
transformed to zero everywhere in 2-D or 3-D models. Ions near z 
= _ 19.5 r e in the nine models used here drift earthward at a speed 
that is about 25% less than E y /B z0 based on B z at z = 0. The speed 
E !B Z0 at x = -19.5 R E is about 140, 120, and 85 km/s in the thick, 
standard, and thin models, respectively. At x = -14.5 /?£, the model 
was only about 10% less than E y /B ZO y which equaled 85, 75, 
and 65 km/s respectively in the three models. We therefore con- 
clude that E y /B zo produces approximate but not accurate estimates 
of V ix in 2-D or 3-D models. 
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Figure 11. The r, y , and z components of the ion bulk velocity for the standard 5-keV model. 


The single particle orbits show the physical origin of the strong 
field-aligned drift. The effect is produced primarily by ions that 
mirror within the current sheet. For example, the ion in Figure 2a 
moved approximately 5 R E earthward, generally along B, as it went 
from x = -19 R E , z = 0.5 R E to a mirror point at * = -14 R E . The ion 
then moved only about 2.5 R E tailward from the mirror point, pri- 
marily along B, until it returned to z = 0.5 R E . The net 2.5 R E earth- 
ward motion at z > 0.5 R E contributes to a generally field-aligned 
bulk flow. Other examples are evident in Figures 2c and 3a. 

A small net drift toward z = 0 is seen in the model V iz . This is 
primarily E x B drift in the negative z direction, though 
field-aligned drift also provides a positive contribution to V z . The 
model E y = 0.3 mV/m electric field produces a 20 km/s drift toward 
z = 0 in a region with B x = 15 nT. Electrons and ions E x B drift in 
the same direction, so these contributions to the current tend to 
cancel in the outer current sheet, where the guiding center approxi- 
mations are valid. 

Although it does not affect the COT analysis, which is based 
only on j yi electrons were assumed to drift along field lines at the 
same speed as ions. This assumption was needed to make 
V * j =0 for electrons, and also so there would be no Birkeland 
current. No field aligned current is expected in the present steady 
state model with essentially no y dependence in the region of inter- 
est. 

6. Discussion and Summary 

The Consistent Orbit Tracing (COT) technique was described. 
The principal steps were: select a number of proton groups, trace 
orbits of the 1000 protons in each group, evaluate the cross-tail 
current carried by each group of ions plus their associated elec- 


trons, and combine these groups so that they produce a nearly con- 
sistent current sheet. In a consistent current sheet, the ions and 
electrons carry the current needed to generate the magnetic field in 
which the orbits were traced. 

A set of nine COT analyses were carried out using 3 proton 
energies (1.5, 5, and 15 keV) for each of three magnetic field mod- 
els. Characteristic current sheet thickness scale lengths varied from 
0.3 R e to 1 R e in the models used. Selection of these thicknesses 
was guided by the few available observations of the instantaneous 
structure of a steady magnetotail current sheet. The adiabaticity 
parameter k ranged from 0.3 to 0.65 in the -20 R E < x < -14 R E 
region of interest. A uniform zeroth order cross-tail electric field 
was present in all calculations. A first order parallel electric field 
was added to maintain charge neutrality. 

An earlier study using a 1-D modified Harris magnetic field 
model [Kaufmann and Lu t 1993] was unable to generate a current 
sheet that agreed with observations in the region of interest. The 
introduction of a 2-D equilibrium tail module, which has both x 
and z dependence, removed the inconsistencies found in the earlier 
1-D analysis. It was possible to generate nearly consistent current 
sheets for most of the magnetic field model and particle energy 
combinations. Adequate solutions were not found when attempting 
to combine the thicker model current sheets and ions with energies 
that were lower than those typically observed in the region of inter- 
est. 

The results of this analysis support the postulate that in 
sheet-like structures, current is carried mostly by particles with 
B zo^xo < K < 1 . where B xo is B x in the lobes and B zo is B z at z = 
0. The lower limit on k is needed for consistency in a modified 
Hams 1-D model field. It is based on the need to confine current 
within the current sheet in any consistent model. In the present 
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study we found it difficult to create a consistent current sheet using 
particles with k > 1 . The ring current magnetic field is generated 
primarily by particles inside the Earth, so k > 1 particles can domi- 
nate in the radiation belts. 

Burkhart et al. [1992] carried out a consistent orbit tracing study 
for a 1-D magnetotail. Their treatment was different from that used 
here in that it involved a fixed source of particles injected at the 
edge of the current sheet. The magnetic field was allowed to 
change until consistency was attained. The model therefore had 
few trapped particles. The resulting current sheets were very thin, 
with small k and a scale size of approximately z 0 . Burkhart et al 
[1992] concluded that this type of current sheet could not be gener- 
ated when k approached 1 . 

Ashour-Abdalla et al. [1993, 1995] used a 2-D version of the 
T89 model for a series of orbit tracing studies. They used a fixed 
source in the mantle for these large-scale kinetic (LSK) models. A 
goal was to demonstrate that a reasonable magnetotail could be 
generated using only a mantle source. The particles were followed 
as they drifted all the way in from the distant tail through our 
region of interest. There were so many chaotic interactions with the 
current sheet that all orbit categories became populated by the time 
particles reached our region of interest. 

The final model was found to be nearly in force balance, even 
though the model was not really consistent because both the parti- 
cle source and the basic fields were preselected. Calculated cur- 
rents tended to be carried in a sheet that is thinner than would be 
needed to generate the T89 field. The density profile was qualita- 
tively similar to our Figure 10a, though the drop-off in z was 
slower than in our results. Earthward drift was slower than in our 
model because the LSK calculations were run with an 0.1 mV/m 
cross-tail field, while our COT calculations used 0.3 mV/m. The 
cross-tail drift velocities were comparable in the two models, as is 
required to support the comparable lobe fields. 

A principal goal of the present study was to use individual parti- 
cle orbits to better understand how particles carry current in a 
sheet-like structure. All ions in each group were started with the 
same energy so that the group would be dominated by a particular 
type of orbit. Particle orbits previously have been classified accord- 
ing to the dynamical properties of the particles. Guiding center, 
chaotic, and resonant orbits are examples. Instead, orbits in this 
work were grouped according to the characteristic spatial distribu- 
tion of cross-tail current carried by the particles. This resulted in 
three orbit categories: particles that remain near z = 0 (trapped 
orbits), particles that spiral and then magnetically mirror elsewhere 
in the principal current sheet (cucumber orbits), and particles that 
magnetically mirror closer to Earth (Speiser orbits). Consistent 
current sheets required the use of all three orbit categories in each 
of the nine cases studied. 

A study of single ion orbits showed the features that resulted in 
the characteristic j y (x, z) pattern earned by particles in each orbit 
category. Magnetization currents produced most of the structure 
and characteristic scale of the magnetotail current sheet. Substan- 
tial earthward drift was found in the COT results. This was attrib- 
uted to particles mirroring within the current sheet. Such particles 
were seen to drift farther earthward than tailward along field lines 
during a bounce period. 

A number of assumptions were introduced in order to include 
electrons. Some of these assumptions were needed to satisfy basic 
requirements of the steady state model, such as charge neutrality 
and conservation. Other assumptions, such as the ratio of electron 
to ion energy, were based on observations. Electrons with the ener- 
gies used here obeyed the guiding center approximations through- 
out the region of interest in all three magnetic field models. 


Cross-tail electron currents were relatively small in the model 
current sheets studied here. Electron current can be much more 
important in a thinner current sheet. Extremely thin current sheets 
have been observed just before and just after substorm onsets 
[Mitchell et al., 1990; Sergeev et al ., 1993; Hesse et al., 1995; 
Pritchett and Coroniti , 1995]. The present work did not attempt to 
model such situations. 

The COT analysis produced an ion distribution function for 
each of the 120 spatial boxes used in this study. The distribution 
functions were integrated to evaluate fluid parameters. Three of the 
simplest parameters: the ion density, temperature, and bulk flow 
velocity, were in good agreement with observations. Comparisons 
of other fluid parameters, a study of force balance, and an analysis 
of energy fluxes will be presented in the future. 
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